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Simulation  of  Oil  Slick  Transport  in 
Great  Lakes  Connecting  Channels 
Theory  and  Model  Formulation 

HUNG  TAO  SHEN,  POOJITHA  D.  YAPA  AND  MARK  E.  PETROSKI 


INTRODUCTION 

In  recent  years  there  has  been  a  growing  concern  over  the  increasing  contamination  of  waterways 
and  adjacent  shoreline  areas  caused  by  oil  spills.  The  Great  Lakes-St.  Lawrence  River  system  (Fig. 
1)  is  the  world's  largest  freshwater  system.  Being  a  navigation  corridor,  it  is  subject  to  the  risk  of  oil 
contamination  (Beurket  and  Argiroff  1984,  Argiroff  and  Weigum  1986).  With  the  possibility  of  oil 
spills  in  the  system,  an  adequate  means  is  needed  for  analyzing  or  predicting  the  movement  and 
spreading  of  oil  slicks.  In  this  study,  computer  simulation  models  were  developed  for  oil  slick  trans¬ 
port  in  lakes  and  rivers.  These  models  can  either  be  used  on  a  real-time  basis  to  predict  the  movement 
of  an  oil  spill  to  assist  the  clean-up  or  as  scenario  models  to  analyze  possible  impact  of  oil  spills  from 
navigation. 


UNITED  STATES 

Figure  1.  Great  Lakes-St.  Lawrence  River  System. 


The  programs  are  written  in  Fortran  programming  language  to  be  compatible  with  a  Fortran77 
compiler.  In  addition,  a  user-friendly,  menu-driven  program  with  graphics  capability  was  devel¬ 
oped  for  the  IBM-PC  AT  computer,  so  that  these  models  can  be  easily  used  to  assist  the  oil  spill  clean¬ 
up  should  a  spill  occur.  This  report  provides  a  complete  description  of  the  analytical  formulation  of 
the  models.  The  logic  and  structures  of  the  computer  programs  and  the  instructions  for  using  the 
models  are  described  in  separate  reports  (Shen  et  al.,  in  prep,  a,  b;Yapa  et  al.,  in  prep.). 

Fate  of  oil  slicks 

To  develop  an  oil  spill  model,  it  is  necessary  to  understand  the  transformation  process  of  an  oil 
slick.  In  addition  to  the  location,  size  and  physical-chemical  properties  of  the  spill,  the  fate  of  an  oil 
slick  in  a  water  body  is  governed  by  complex  interrelated  transport  and  weathering  processes.  Table 
1  summarizes  the  environmental  pathways  of  a 
typical  crude  oil  spill  at  sea.  Immediately  upon 
entering  a  water  body,  the  spilled  oil  spreads 
and  forms  a  surface  slick  that  covers  a  large  area 
of  the  water  surface  and  can  be  moved  about  by 
the  action  of  winds,  waves  or  currents.  Some 
light  hydrocarbons  and  some  polar  components 
begin  to  go  into  solution  in  the  underlying  water 
column,  but  most  of  these  are  lost  to  the  atmo¬ 
sphere  through  evaporation.  Evaporation  of  vola¬ 
tile  components  and  dissolution  of  hydrocar¬ 
bons  may  reduce  the  volume  of  spilled  oil  by  as 
much  as  50%  during  the  first  few  days  following 
the  spill.  In  turbulent  waters,  some  of  the  oil  is 
emulsified  into  the  water  column  as  small  dispersed  droplets.  These  droplets  may  become  dispersed 
by  currents,  or  they  may  become  attached  to  suspended  particulate  matter  and  slowly  settle  to  the 
bottom.  The  turbulence  action  can  also  cause  water  to  become  entrained  in  the  oil,  forming  water-m¬ 
oil  emulsion  that  may  eventually  weather  further  to  form  dense  tar  balls.  While  all  these  processes 
are  occurring,  photochemical  reactions  and  microbial  biodegradation  can  also  change  the  character 
of  the  oil  and  reduce  the  amount  of  oil  present. 

These  physical  and  chemical  processes  can  be  categorized  as  advection,  spreading,  evaporation, 
dissolution,  emulsification,  auto-oxidation,  sedimentation  and  biodegradation.  A  number  of  excel¬ 
lent  reviews  on  these  processes  have  been  published  (Stolzenbach  et  al.  1977,  Wheeler  1978,  Jordan 
and  Payne  1980,  Huang  and  Monastero  1982).  Figure  2  is  a  schematic  representation  of  the  transport 
and  weathering  processes  showing  time  periods  for  which  each  of  these  processes  is  important.  It  is 
clear  that  the  oil  undergoes  significant  modification  during  the  first  few  hours  following  a  spill.  A 
brief  discussion  of  each  of  these  physical,  chemical  and  biological  processes  follows. 

Adwction 

Advection  is  a  physical  process  that  involves  the  drifting  of  the  surface  oil  slick  and  the  subsurface 
oil.  The  advection  of  surface  oil  is  caused  by  the  combined  effects  of  surface  current  and  wind  drag. 
The  advection  of  subsurface  oil  is  the  movement  of  the  oil  entrained  in  the  flow  due  to  the  subsurface 
current. 

Spreading 

Spreading  also  involves  the  areal  expansion  of  oil  slicks.  The  spreading  process  includes  mechani¬ 
cal  spreading  and  turbulent  diffusion.  In  mechanical  spreading  the  area  of  the  surface  oil  slick  ex¬ 
pands  because  of  the  balancing  forces  of  inertia,  gravity,  viscosity  and  surface  tension.  T urbulent  dif¬ 
fusion  is  the  spreading  of  the  surface  oil  because  of  the  fluctuations  of  the  wind  and  current  velocities. 


Table  1.  Pathways  for  the  environmc  'al  fate 
of  crude  oil  (Butler  et  al.  1976). 

Time  scale  Percentage  of 
Pathway  Ulays)  initial  oil 


Evaporation 

1-10 

25 
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1-10 
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15 

Residue 

>100 
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100 
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Figure  2.  Physical,  chemical  and  biological  processes  affecting  oil  slick  transformation. 


Since  the  spreading  of  oil  will  increase  weathering  processes  such  as  evaporation,  dissolution,  photo¬ 
oxidation  and  biodegradation,  it  is  one  of  the  most  important  processes  affecting  the  fate  of  the  spilled 
oil. 

Evaporation 

Evaporation  begins  immediately  after  the  spill.  As  the  spill  spreads,  more  of  the  hydrocarbons  are 
exposed  to  the  atmosphere,  causing  evaporation  rates  to  increase.  The  amount  and  rate  of  evapora¬ 
tion  depend  on  the  percentage  of  light  (or  volatile)  components  in  oil.  Evaporation  is  the  most  sig¬ 
nificant  physical-chemical  process  in  reducing  the  oil  volume.  Highly  refined  oils  can  lose  75%  or 
more  of  their  volume  through  evaporation  within  days. 

Dissolution 

Some  oil  components  will  be  lost  from  a  slick  by  dissolving  in  the  water  column.  Only  the  low- 
molecular-weight  hydrocarbons  have  an  appreciable  solubility  in  water.  The  fraction  of  oil  dissolved 
is  very  small  compared  to  that  evaporated,  but  this  extraction  process  can  be  important  because  of 
the  toxicity  of  the  dissolved  fraction. 

Emulsification 

The  formation  of  oil-in-water  emulsion  is  an  important  mode  of  oil  dispersion.  Turbulence  and 
wave  action  mix  the  surface  layer  of  oil  into  the  water  column  by  causing  many  very  small  globules 
of  oil  to  form,  which  can  be  rapidly  dispersed  vertically  into  the  water  column  where  they  are  subject 
to  subsurface  transport.  The  eventual  fate  of  oil-in-water  emulsions  is  probably  dissolution  in  the 
water  column,  attachment  to  solid  particles  and  biodegradation.  Water-in-oil  emulsions  can  also 
form,  particularly  with  heavy  crudes  and  residual  oils.  The  resulting  emulsions,  which  contain  a 
large  percentage  of  water  but  have  a  semisolid  texture,  are  often  referred  to  as  "chocolate  mousses" 
because  of  their  appearance.  These  emulsions  may  persist  on  the  water  surface  and  disintegrate  into 
tar  lumps  after  a  long  time  if  they  do  not  wash  up  on  the  shore. 
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Auto-oxidation 

In  the  presence  of  atmospheric  oxygen,  natural  sunlight  has  sufficient  energy  to  change  the  com¬ 
position  of  the  oil.  This  photo-oxidation  process  is  very  slow.  The  extent  and  rate  of  the  photo-oxida¬ 
tion  primarily  depend  on  the  chemical  composition  of  the  oil  and  its  optical  density.  Little  is  known 
about  the  effect  of  photochemical  reactions  in  the  overall  oil  slick  transformation  process. 

Sinking  and  sedimentation 

Some  of  the  spilled  oil  ultimately  sinks  if  it  does  not  wash  up  on  the  shore.  This  sinking  and  sedi¬ 
mentation  process  occurs  because  the  density  of  the  oil  increases  when  the  lighter  fractions  of  the  oil 
evaporate  or  dissolve  or  when  the  oil  adheres  to  suspended  sediment.  This  process  may  eventually 
cause  oil  fractions  to  sink  to  the  bottom,  where  they  may  be  moved  laterally,  they  may  be  resus¬ 
pended,  or  they  may  undergo  further  biological  or  physical-chemical  reaction. 

Little  is  known  about  the  ultimate  fate  of  the  sedimented  oil.  All  of  the  processes  just  described, 
except  possibly  photo-oxidation,  can  only  redistribute  the  oil.  They  cannot  remove  the  hydrocarbon 
from  the  environment.  Real  degradation  takes  place  through  biochemical  oxidation;  this  process  is 
the  principal  long-term  means  of  removing  the  spilled  oil  from  the  environment. 

Oil  spill  simulation  models 

Many  of  the  oil  spill  models  developed  during  the  last  decade  simulate  only  the  advection  and 
spreading  processes.  Other  models  deed  extensively  with  physical-chemical  processes  but  lack  the 
component  for  simulating  the  advection  of  the  slick.  Only  in  recent  models  have  the  incorporation 
of  both  the  transport  and  weathering  processes  been  attempted  (Huang  andMonastero  1982).  Since 
there  are  not  enough  data  to  establish  a  reliable  analytical  formulation  for  many  of  the  weathering 
processes,  it  is  impractical  to  include  all  of  them  in  an  oil  spill  simulation  model.  It  would  be  more 
useful  to  include  the  most  significant  processes  (i.e.,  those  accounting  for  the  bulk  of  the  oil)  while 
omitting  others  so  that  uncertainty  in  the  outcome  can  be  reduced.  In  addition,  since  part  of  the  oil 
spilled  in  water  bodies,  especially  in  rivers  and  lakes,  will  be  washed  up  on  shore,  appropriate 
shoreline  boundary  conditions  must  also  be  considered  in  a  simulation  model. 

Almost  all  of  the  existing  oil  slick  models  were  developed  for  coastal  marine  environments.  Only 
a  few  models  were  developed  for  rivers  and  lakes  (Huang  and  Monastero  1982).  Tsahalis  (1979b)  de¬ 
veloped  a  simulation  model  for  predicting  the  transport,  spreading  and  associated  shoreline  contam¬ 
ination  of  oil  spills  in  rivers.  In  this  model  the  current  velocity  distribution  in  the  river  is  calculated 
by  empirical  relationships  determined  from  field  data,  with  some  modifications  for  the  secondary 
current  in  riverbends.  The  oil  slick  is  assumed  to  remain  circular,  with  its  radius  calculated  according 
to  Fay* s  spreading  laws  (Fay  1 969, 1 971 ).  The  drift  veloci  ty  of  the  slick  is  determined  by  formulas  de¬ 
rived  by  Tsahalis  (1979)  from  laboratory  experiments. 

Fingas  and  Sydor  (1980)  developed  a  two-dimensional  model  for  oil  spills  in  rivers.  In  this  model 
the  current  velocity  distribution  is  determined  by  the  two-dimensional  finite-difference  scheme  of 
Leendretse  (1970) .  The  entire  oil  slick  volume  is  represen  ted  by  a  large  number  of  individual  particles. 
The  drift  velocities  of  these  particles  are  determined  by  the  wind  factor  approach.  A  random  fluc¬ 
tuation  component  is  included  to  represent  horizontal  diffusion.  The  spreading  of  the  oil  slick  is  cal¬ 
culated  by  Fay's  spreading  laws  for  circular  slicks.  For  oil  spills  in  lakes,  the  only  model  is  one  de¬ 
veloped  at  the  Great  Lakes  Environmental  Research  Laboratory  (Boyd  1979,  Schwab  et  al.  1984)  for 
the  Great  Lakes;  this  is  basically  a  model  for  predicting  the  motion  of  a  group  of  surface  particles.  The 
movement  of  an  oil  slick  can  be  simulated  using  this  model  by  representing  the  oil  slick  as  a  group 
of  particles.  None  of  these  three  models  considered  oil  weathering  processes,  so  the  effects  of 
weathering  on  the  oil  slick  transformation  cannot  be  accounted  for. 

In  this  study,  computer  models  were  developed  for  simulating  the  fate  of  oil  spills  in  a  river  or  lake 
including  the  effect  of  ice  covers.  None  of  the  previous  models  for  oil  spills  in  rivers  and  lakes  took 
this  factor  into  consideration.  The  purpose  of  these  simulation  models  is  to  assist  the  on-scene  com- 
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INPUT  DATA 

TYPE  OF  SPILL 

VOLUME  AND 

TYPE  OF  OIL 

LOCATION 

RIVER/LAKE  DATA: 

1)  Instantaneous 

Density,  Viscosity, 

AND  TIME 

OF  SPILL 

Wind  Velocity 

2)  Continuous 

Surface  Tension,  etc. 

Discharge  and  Water  Levels 

Temperature 

Ice  Condition 

DETERMINE  WATER  CURRENT 
Two-dimensional  Depth-averaged 
Velocity  Distribution 


ADVECTION 


u 


MECHANICAL  SPREADING 
AND 

TURBULENT  DIFFUSION 


SHORELINE  DEPOSITION 


WEATHERING 
(Evaporation,  Dissolution) 


OUTPUT  DATA 


’  Two-dimensional  Distribution  of  Oil  in  the  River 
'  Amount  and  Locations  of  Deposition  on  Shoreline 


Figure  3.  Structure  of  the  simulation  model. 

mander  to  develop  clean-up  measures  in  the  case  of  an  actual  spill  and  to  assess  likely  environmental 
impacts  of  possible  spills.  The  models  are  primarily  designed  for  predicting  the  motion  of  an  oil  slick 
in  rivers  or  lakes. 

Figure  3  shows  the  structure  of  the  simulation  model.  Detailed  discussions  of  the  model  formu¬ 
lation  will  be  presented  later.  The  models  developed  in  this  investigation  consider  the  oil  slick  to  be 
composed  of  a  large  number  of  discrete  parcels,  which  are  tracked  for  their  positions  and  volume  to 
each  time  level  during  the  period  of  simulation.  The  two-dimensional  velocity  distributions  of  the 
underlying  river  or  lake  water  are  first  computed.  Using  the  wind  velocity  and  the  computed  current 
velocity,  the  model  determines  the  advection  of  each  parcel  in  the  slick  using  the  wind  factor  ap¬ 
proach.  The  spreading  of  the  slick  is  simulated  by  considering  both  the  mechanical  spreading  and 
surface  turbulent  diffusion.  The  model  developed  by  Fay  (1969, 1 971 )  is  used  to  simulate  mechanical 
spreading,  while  a  random-walk  simulation  is  used  to  account  for  the  spreading  due  to  surface  turbu- 
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lence.  The  loss  of  oil  due  to  shoreline  deposition  is  calculated  according  to  the  oil  retention  capability 
of  the  shoreline  that  the  oil  slick  reaches.  The  loss  of  oil  due  to  evaporation  and  dissolution  are  cal¬ 
culated  based  on  empirical  formulations  that  consider  the  effects  of  slick  area,  wind  velocity, 
temperature  and  oil  properties.  Since  weathering  processes  that  occur  long  after  the  onset  of  spill  are 
not  well  understood  and  are  less  significant,  these  processes  are  not  considered  in  the  model.  This  is 
justified  from  the  operational  point  of  view,  since  the  oil  will  be  washed  up  on  shore  and  clean-up 
measures  will  beg'.i  shortly  after  the  spill. 


THEORY  AND  MODEL  FORMULATION 


Analytical  framework 

The  transformation  of  an  oil  slick  is  affected  by  a  number  of  complex  physical  and  chemical  pro¬ 
cesses.  To  facilitate  the  discussion  of  the  model  formulation,  an  analytical  framework  based  on  the 
equations  of  motion  of  a  surface  oil  slick  will  first  be  presented.  For  a  surface  oil  slick,  as  shown  in 
Figure  4,  the  two-dimensional  depth-averaged  equation  of  motion  can  be  written  as  (Ahlstrom  1975, 
Stolzenbach  et  al.  1977) 


d(pjl)  Xpouh) )  d(Povh)  „ 

dt  dx  dy  Vs  vb 

(1) 

3u+„aw+l,au-_JPw-PoP  _*bx+*« 
dt  dx  dy  °  \  p  /  dx  p/i  ph 

(2) 

dv  +udv  +vdv  =  (Pw-Po)ah  _!by+!^ 

dt  dx  dy  \  P  '  dy  ph  ph 

(3) 

where  x  and  y  =  space  variables 

t  =  time 

u  and  v  =  components  of  the  oil  velocity 
po  and  pw  =  densities  of  oil  and  water,  respectively 


h 

Z 


x.  and  x, 

bx  by 

x  and  x 

sx  sy 


*b 

R 


oil  slick  thickness 
acceleration  due  to  gravity 

shear  stress  components  on  the  bottom  of  the  oil  slick 
shear  stress  components  on  the  top  surface  of  the  oil  slick 
rate  of  oil  flux  outward  through  the  top  surface,  e.g.,  evaporation 
rate  of  oil  flux  outward  through  the  bottom  surface,  e.g.,  dissolution 
other  source  and  sink  terms. 


Wind 


6 


The  density  of  oil  can  be  determined  from  its  specific  gravity.  Typical  values  of  the  specific  gravity 
are  0.7  for  gasoline,  0.8  for  kerosene,  0.84  for  diesel  or  no.  2  fuel  oil,  and  0.98  for  Bunker  C  or  no.  6  fuel. 
The  specific  gravity  of  crude  oils  varies  between  0.80  and  0.99  (Bishop  1983). 

The  boundary  condition  at  the  circumference  of  the  slick  is  that  the  force  parallel  to  the  water  sur¬ 
face  and  normal  to  the  slick  boundary  f  be  balanced  by  the  net  surface  tension  force: 

f  =  c  =  <y  -a  cosG  -  a  cosQ 

J  n  aw  oa  on  ow  ow 


where  a  =  surface  tension  of  the  air/water  interface  (72.75  dynes/cm  at  20 °C) 
om  =  surface  tension  of  the  air/oil  surface  (20  dynes/cm  for  many  crudes) 
a  =  surface  tension  of  the  oil/water  interface  (varies  between  15  and  25  dynes/ cm) 

®oa ,  0)h,  =  contact  angles  between  oil  and  air  and  between  oil  and  water  (Stolzenbach  et  al.  1977). 

Equations  1-4  clearly  show  that  the  movement  of  an  oil  slick  is  governed  by  the  advection  due  to  the 
viscous  forces  acting  on  the  top  and  bottom  surfaces  of  the  slick,  the  spreading  due  to  gravitational, 
viscous  and  surface  tension  forces,  and  the  changes  in  the  mass  and  physical-chemical  properties  of 
the  oil  due  to  various  weathering  processes. 

A  number  of  studies  have  attempted  to  analyze  in  detail  the  hydrodynamic  problem  defined 
above  (Kerr  and  Babu  1970,  DePietro  and  Cox  1979,  Foda  and  Cox  1980).  These  complex  mathemati¬ 
cal  treatments  are  not  suitable  for  real  field  problems.  In  the  present  study,  a  Lagrangian  discrete- 
parcel  algorithm  is  used.  In  this  algorithm  the  oil  slick  is  viewed  as  a  large  ensemble  of  small  parcels. 
Each  parcel  has  associated  with  it  a  set  of  spatial  coordinates  and  a  specific  quantity  of  mass.  The 
movement  of  each  parcel  in  the  river  or  the  lake  is  affected  by  the  wind,  water  currents  and  the  con¬ 
centration  of  surrounding  parcels.  If  a  large  number  of  parcels  are  released  in  a  water  body  and  their 
discrete  paths  and  masses  are  followed  and  recorded  as  functions  of  time  relative  to  a  fixed  reference 
grid,  then  the  density  distribution  of  the  ensemble  in  the  water  body  can  be  interpreted  as  the  con¬ 
centration  of  the  oil.  This  approach  requires  an  efficient  bookkeeping  procedure  rather  than  the  so¬ 
lution  of  a  large  matrix  associated  with  a  conventional  Eulerian  finite-difference  or  finite-element 
methods.  The  algorithm  is  inherently  stable  with  respect  to  time  steps,  although  the  time  step  should 
be  compatible  with  the  grid  size  and  velocity  for  numerical  accuracy  (Cheng  et  al.  1984).  Since  the 
movement  of  each  parcel  in  the  oil  slick  depends  on  the  distribution  of  the  entire  ensemble,  all  parcels 
must  be  traced  at  each  time  level  before  proceeding  to  the  next. 

The  detailed  structure  and  implementation  of  the  numerical  model  will  be  discussed  later.  In  the 
following  sections  the  analytical  formulations  for  each  component  of  the  model  will  be  discussed. 

River  current  simulation 

Since  the  water  current  affects  both  the  advection  and  the  spreading  of  an  oil  slick,  the  distribution 
of  both  the  magnitude  and  the  direction  of  the  current  must  be  determined  first.  There  are  numerous 
numerical  methods  for  determining  the  two-dimensional  flow  distribution  in  shallow  waters 
(Leendretse  1970,  Grubert  1976,  Hamilton  et  al.  1982).  However,  because  of  stability  and  accuracy 
problems  of  the  numerical  method,  the  geometric  element  size  A  and  the  size  of  the  time  steps  t  are 
limited  by  t  <,  {AJghf /2.  This  makes  the  method  impractical  for  long  rivers.  In  view  of  the  need  for 
simulating  flow  distributions  in  long  river  reaches,  an  alternative  approach  is  taken  in  the  present 
study.  In  the  present  approach,  the  time-dependent  discharge  distribution  Q(x,t)  along  the  river  is 
first  obtained  by  using  a  one-dimensional  hydraulic  transient  model,*  which  was  developed  based 
on  the  St.  Venant  equations: 

•For  example,  the  unsteady  flow  model  for  the  Great  Lakes  connecting  channels  developed  by  R.  Thomas  of  the  Detroit 
District,  US.  Army  Corps  of  Engineers. 
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*2+^i  =  o 

3x  3 t 


(5) 


dQ 
3 t 


dx  A  dx  ox 


where  x 


A 

H 

S 

o 

Sf 


=  longitudinal  distance  along  the  river 
=  flow  cross-sectional  area 
=  water  level 
=  bed  slope 
=  frictional  slope. 


(6) 


The  friction  slope  can  be  calculated  by  Manning's  equation: 


Sf  = 


" h^Q2 

2.21  A2R4'3 


(7) 


where  R  is  the  hydraulic  radius  and  «b  is  the  Manning's  roughness  coefficient  for  the  bed.  For  an  ice- 
covered  reach,  the  composite  Manning's  coefficient,  which  accounts  for  the  resistance  of  ice  cover  and 
the  riverbed,  should  be  used  instead  of  n  fc.  The  composite  Manning's  coefficient  n  can  be  calculated 
by  the  Belokon-Sabaneev  formula: 

«  =  [l(n,3/*wg/2  )]2/3  (8) 

where  n.  is  the  Manning's  roughness  coefficient  for  the  ice  cover.  Tire  hydraulic  radius  can  be  as¬ 
sumed  to  be  half  of  the  flow  depth  for  ice-covered  reaches. 

Once  the  one-dimensional  solution  is  obtained,  the  discharge  Q  can  be  distributed  across  the 
width  of  the  river  using  a  stream-tube  model  (Shen  and  Ackermann  1980)  to  give  the  two-dimension¬ 
al  velocity  distribution  at  selected  cross  sections. 


Two-dimensional  velocity  distribution 

For  any  channel  cross  section,  as  shown  in  Figure  5,  the  transverse  distribution  of  the  flow  can  be 
determined  using  a  simplified  method  developed  by  Shen  and  Ackermann  (1980).  In  this  method, 
the  cross  section  is  first  divided  into  trapezoidal  elements.  By  applying  Manning's  equation  to  the 
ratio  of  discharges  between  the  entire  cross  section  and  a  partial  cross  section  (Fig.  5),  the  following 


a.  Cumulative  discharges  over  a  b.  Area  and  discharge  for  the  p"'  trapezoid, 

partial  cross  section. 

Figure  5.  Method  for  determining  the  transverse  distribution  of  the  flow  for  any  cross  section. 
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expression  can  be  written: 


n  2 
Qp  _  p=i 


Vp 


2/3 


N 

I  AnRn 

ii=i 


2/3 


where  P  =  number  of  trapezoids  in  the  partial  cross  section 

N  =  total  number  of  trapezoids  describing  a  cross  section  geometry 
Qp  =  cumulative  discharge  up  to  and  including  the  p‘h  trapezoid 
Q  =  the  total  discharge  through  the  entire  cross  section 
An  and  Rn  =  area  and  hydraulic  radius,  respectively,  of  the  «lh  trapezoid 
Ap  and  Rp  =  area  and  hydraulic  radius,  respectively,  of  the  plh  trapezoid. 


The  cumulative  discharge  Qp  can  be  computed  by  first  rewriting  eq  9  as 

Qp  =  Fq  l  [ApRp2/y) 

P=l 

where 


Fq  = 


Q 


2  (^n273) 

n= 1 


and 


(9) 


(10) 


(11) 


Qp-Qp_i-FQ-ApRp2/?  (12) 

Based  on  the  computed  distribution  of  Q,„  stream-tube  boundaries  within  the  cross  section  can  be 
determined  by  simple  interpolation.  Once  the  stream-tube  boundaries  are  located,  the  flow  through 
each  stream  tube  Qsis  divided  by  the  cross-sectional  area  of  the  stream  tube  As  to  obtain  the  depth- 
averaged  velocity: 

Vp  =  QsMs.  (13) 

This  velocity  is  then  assigned  to  the  center  of  the  stream  tube.  Applying  this  procedure  to  successive 
cross  sections  along  the  river  produces  a  two-dimensional  depth-averaged  velocity  distribution.  The 
directions  of  velocity  vectors  are  the  same  as  the  vectors  connecting  center  points  of  each  stream  tube 
at  successive  cross  sections.  As  an  example,  the  simulated  depth-averaged  velocity  distribution  for 
a  reach  of  the  St.  Clair  River  is  shown  in  Figure  6a. 

Once  the  velocity  distribution  along  stream  tubes  is  established,  the  distribution  of  velocity  for  all 
the  points  in  a  predefined  grid  system,  as  shown  in  Figures  6b  and  c,  can  be  obtained  through  linear 
interpolations.  Calculated  velocities  compare  favorably  with  data  obtained  in  the  field  (U.S.  Army 
Corps  of  Engineers,  Detroit  District  1974). 

Lake  circulation 

The  following  equations  of  motion  describe  the  circulation  in  shallow  lakes,  if  nonlinear  convec¬ 
tive  terms  are  neglected: 

=sD^L  +  ±-  (Tl-xb-Ti)  (14) 

at  dx  Pw 
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b.  Superimposed  grid  system. 


Figure  6.  Method  for  determining  the  velocity  distribution  fora  typical  river 

and 


dH  +dM  +dN  _Q 
dt  dx  dy 

where  M  =  uD 
N  =  v(D 

uc  and  vc  =  depth-average  velocity  components 
D  =  depth  of  the  flow 


H  =  water  surface  elevation 
g  =  acceleration  due  to  gravity 
/  =  Q  sin<J>  (the  coriolis  parameter) 

£2  =  angular  velocity  of  the  earth 
=  latitude 

and  r®  =  components  of  the  wind  stress 
x£  and  x^  =  components  of  the  bed  shear 

xjj  and  x^,  =  components  of  the  shear  stress  at  the  ice/water  interface 
randy  =  space  variables 
t  =  time. 


For  simulating  the  circulation  pattern  in  the  Lake  St.  Clair,  the  RLID  finite-difference  model  devel¬ 
oped  by  Schwab  et  al.  (1981 )  is  used.  This  model  neglects  the  free  surface  fluctuation.  Equation  1 6  then 
leads  to  the  existence  of  the  stream  function  vy,  defined  by 


M  =  -  N  =  .  (17) 

3y 

Combining  this  with  eq  14  and  15  will  yield 

1  (y.D~1Vv)+/(a¥  dP"1. -dv  dP."1  \-  3  (Ty~Ty~Ty]_  d  (Tx~Tx~Tx)  08) 

dt  \3r  9y  3 y  dx  )  dx\  Dpw  )  dy  \  Dpw  / 


Equation  1 8  is  used  to  solve  for  \p  values  through  a  second-order  finite-difference  scheme.  The  veloc¬ 
ity  distribution  is  determined  from  the  calculated  values.  Figure  7  shows  sample  results  for  the 
velocity  distribution  in  Lake  St.  Clair. 


Figure  7.  Velocity  distribution  in  Lake  St.  Clair.  The  stage  and  discharge  at  the  outlet  of  the  lake  are  575.04  ft 
and  210  POO  ft*/s,  respectively. 


11 


Advection 

Advection  in  open  water 

Advection  or  drifting  of  oil  on  the  water  surface  is  the  result  of  the  combined  action  of  wind,  cur¬ 
rent  and  wind-generated  waves  (Schwartzberg  1971,  Tsahalis  1979a).  The  drift  velocity  of  the  surface 
oil  is  usually  considered  to  be  the  vector  sum  of  a  wind-induced  drift  and  a  water-current-inJuced 
drift  (Stolzenbach  et  al.  1977).  In  the  present  model  the  advective  velocity  of  each  oil  parcel  is 
computed  as 

Vt  =  V  +  V'  (19) 

where  Vt  is  the  drift  velocity  of  an  oil  parcel  and  V  and  V  are  the  mean  and  turbulent  fluctuation 
components  of  the  drift  velocity.  The  component  V'is  included  to  simulate  the  horizontal  diffusion 
of  the  oil  parcels.  The  formulation  for  V'  will  be  discussed  later  in  this  section.  The  mean  velocity 
component  V  includes  both  the  wind  and  water  current  effects  and  can  be  calculated  by  the  formula 

V=awVw  +acVc  (20) 

where  Kw  =  wind  velocity  at  10  m  above  the  water  surface 
Vc  =  «c  +VC  =  depth-averaged  current  velocity 
a  =  factor  to  account  for  the  drift  of  the  surface  slick  due  to  wind 

W 

a,  =  factor  to  account  for  the  drift  of  the  oil  slick  at  the  water  surface  due  to  the  current. 


Wind-induced  surface  currents  have  been  reported  to  vary  between  1  %  and  6%  of  the  wind  speed, 
with  3%  being  the  most  widely  used  drift  factor  in  oil  slick  trajectory  models  (Stolzenbach  et  al.  1977). 
Based  on  Madsen's  (1977)  Ekman  layer  model,  a  drift  factor  of  0.03  can  also  be  obtained  with 
appropriate  values  of  equivalent  roughness  height  of  the  free  surface  (Huang  and  Monastero  1982). 

If  the  velocity  profile  is  assumed  to  be  logarithmic,  the  surface  velocity  of  the  water  current  can 
be  related  to  the  depth-averaged  current  velocity  by  the  relationship 


£.1  + 

Vc 


u» 

kVc 


(21) 


where  Vs  =  surface  velocity 
u,  =  shear  velocity 
K  =  Karman  constant  (0.4). 


Calculations  based  on  conditions  in  the  connecting  channels  show  that  the  ratio  Vs/Vc  varies  between 
1 .1  and  1 .2,  with  most  of  the  values  equal  to  1 .1 .  In  the  present  simulation,  values  of  0.03  and  1 .1  are 
used  for  aw  and  respectively. 

Advection  under  an  ice  cover 

Oil  spills  under  ice  have  received  little  theoretical  or  laboratory  treatment.  As  suggested  by  the 
experimental  study  of  Cox  and  Schultz  (1981),  an  ice  cover  may  be  classified  into  three  categories 
when  determining  slick  advection:  the  ice  cover  may  be  smooth,  it  may  contain  small  roughness 
elements,  or  it  may  contain  large  roughness  elements.  The  following  discussion  summarizes  the 
results  of  limited  experimental  studies  and  serves  as  the  basis  for  modeling  the  advection  of  oil  in  the 
present  simulation.  i 

Under  smooth  ice  covers  with  no  current,  the  oil  will  rest  at  an  equilibrium  thickness  that  is 
described  by  the  empirical  equation 

\ 

5eq  =  1.67-8.5(Apw)  (22) 

where  8^  =  static  equilibrium  slick  thickness  (cm) 
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A  =  (pw  -  pc)/pw  =  relative  density  difference  between  oil  and  water 
pw  and  po  =  water  and  oil  densities,  respectively  (g/cm). 

An  ice  cover  is  considered  to  be  smooth  when  the  height  of  the  ice  roughness  is  smaller  than  the  equi¬ 
librium  thickness  of  the  oil  8  . 

eq 

The  depth-averaged  current  velocity  at  which  the  oil  just  begins  to  move  along  an  ice  cover  is 
called  the  threshold  velocity  l/(h.  For  a  smooth  cover  the  value  of  U ’ft  was  empirically  determined  to 
be  a  function  of  the  oil  viscosity  poand  is  given  as 


Ulh  =  305.79/(88.68 -po) 


(23) 


where  ti(h  is  in  cm/ sec  and  po  is  in  g/cm-s.  Viscosities  for  crude  and  fuel  oils  fall  in  the  range  of  5-50 
centipoises.*  Typical  values  of  oil  viscosity  can  be  found  in  fluid  mechanics  books  (e.g.  Rouse  1946). 

A  rough  ice  cover  can  retain  the  oil  between  the  roughness  elements.  As  the  current  velocity 
increases,  the  oil  creeps  along  the  upstream  face  of  the  roughness  element  until  it  spills  over  the 
element  and  moves  downstream.  The  threshold  current  velocity  at  which  the  oil  will  move 
downstream  under  a  rough  ice  cover  is  called  the  failure  velocity  U„ : 

un  =  (24) 

where  Oq/w  is  the  oil/water  interfacial  tension.  The  failure  velocity  Ufl  is  the  velocity  above  which 
no  oil  can  be  contained  upstream  of  a  large  roughness  element.  If  the  depth-averaged  current  velocity 
is  less  than  the  threshold  velocity,  the  slick  will  not  advect.  If  the  depth-averaged  velocity  is  greater 
than  the  threshold  velocity,  Uth  or  Ufl,  the  relationship  between  the  current  velocity  and  the  slick 
velocity  is 


with 


K 

0.115Fg  + 1.105 


where  V  =  mean  drift  velocity 
V.  =  current  speed 
K  =  friction  amplification  factor 
Fg  =  slick  densimetric  Froude  number 
g  =  acceleration  due  to  gravity. 


(25) 

(26) 


The  value  of  K  is  a  function  of  the  roughness  height  of  the  cover.  With  the  limited  data  available,  K 
is  assumed  to  vary  linearly  between  1 .0  for  a  smooth  cover  and  2.6  for  an  ice  cover  with  a  Manning's 
roughness  coefficient  n,  of  0.055. 

Horizontal  diffusion 

The  term  V'  in  eq  19  is  included  to  account  for  horizontal  diffusion  caused  by  the  turbulent  fluc¬ 
tuation  of  the  drift  velocity.  For  isotropic  turbulent  diffusion,  the  diffusion  coefficient  Ej  can  be  re¬ 
lated  to  the  magnitude  of  V'  by  the  random  walk  analysis  (Fischer  et  al.  1979),  where 


*lcp»  1.0x10s  g/cm-s  *  2.4  Ib/ft-hr. 
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r 


Y  =  (4£T/5f),/2  (27) 

where  8f  is  a  time  step. 

Murray  (1972),  using  the  Fickian  diffusion  theory,  estimated  the  value  of  the  diffusion  coefficient 
for  the  1970  Chevron  spill  to  be  19  m2/s.  Cole  et  al.  (1973)  used  a  value  of  4.5  ft2/ s,  based  on  dye  tests 
for  the  simulation  of  an  oil  spill  at  the  Cherry  Point  Refinery  in  the  Strait  of  Georgia.  Hunter  (1980), 
by  fitting  observed  spreading  of  oil  slicks  to  Okubo's  (1962)  theory,  suggested  a  value  of  a  5  m2/ s.  All 
of  these  analyses  were  based  on  coastal  oil  spills.  In  rivers  and  shallow  lakes  the  diffusivity  will  be 
affected  by  the  shear  velocity  u.  and  the  depth  of  flow  D,  in  addition  to  the  wind.  Flume  experiments 
by  Sayre  and  Chang  (1969)  indicated  that  £T  is  on  the  order  of  0.6  Du..  This  expression  may  be  used 
to  calculate  turbulent  diffusion  in  rivers.  For  lakes  the  following  formula  for  surface  dye  patches 
(Okubo  1962)  may  be  used: 

£t  =  0.0027  tlM  (28) 

where  ET  is  in  cm2/s  and  f  is  in  seconds. 

In  the  present  simulation  the  fluctuation  velocity  component  V '  is  calculated  by 


V'  =V'  rne'0'.  (29) 

For  a  selected  ET  value  the  magnitude  of  Y  is  computed  using  eq  27.  The  variable  rn  is  a  normally 
distributed  random  number  with  a  mean  value  of  0  and  a  standard  deviation  of  1.  The  directional 
angle  0'  is  assumed  to  be  a  uniformly  distributed  random  angle  ranging  between  0  and  n. 

During  each  time  step  Af,Jhe  displacement  AS  of  each  oil  parcel  is  calculated  by  numerically 
integrating  the  drift  velocity  V t  over  the  time  period  At.  For  numerical  accuracy,  when  the  At  value 
is  large,  subintervals  8f  k  are  used  for  the  advection  of  oil  parcels.  In  this  case  the  displacement  during 
the  time  interval  Af  is 

AS  =  £  Vk5fk  (30) 


where  Vk  =  drift  velocity  of  an  oil  parcel  during  the  time  interval  Sfk 
k  AS  -  displacement  during  the  time  interval  At 
Z5fk  =  At. 

The  values  of  8fk  should  satisfy  the  condition  (Roache  1972,  Cheng  et  al.  1984) 

U*  Ay) 

where  uk  and  n  are  the  x  and  y  components  of  the  velocity  V.. 


Mechanical  spreading 

Spreading  in  open  water 

Spreading  of  the  oil  slick  is  one  of  the  most  important  processes  in  the  early  stage  of  the  oil  slick 
transformation,  because  of  the  influence  of  the  surface  area  of  the  oil  slick  on  weathering  processes 
such  as  evaporation  and  dissolution.  Besides  the  horizontal  dispersion  due  to  advection  and  turbu¬ 
lent  diffusion,  the  spreading  of  an  oil  slick  is  determined  by  the  balance  of  gravitational,  viscous  and 
surface  tension  forces.  Spreading  is  also  affected  by  weathering  processes,  which  tend  to  change  the 
mass  and  physical-chemical  properties  of  the  slick.  Several  models  have  been  proposed  for  the  pro¬ 
cess  of  mechanical  spreading  in  open  water  (Blokker  1964,  Fay  1969, 1971,  Hoult  1972,  Mackay  et  al. 
1980).  In  this  study,  Fay's  spreading  theory  (1971)  is  used  because  this  theory  is  based  on  a  compre¬ 
hensive  description  of  the  spreading  mechanisms  and  has  been  verified  by  laboratory  experiments 
(Hoult  and  Suchon  1970,  Fay  1971)  and  other  analytical  solutions  (Fannelop  andWaldman  1972). 


L 
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Table  2.  Spreading  laws  for  oil  slicks  (Fay  1971,  Hoult  1972, 
Waldman  et  al.  1973). 


Spreading  phase  L,  (one-dimensional)  R(  ( radial) 

Gravity-inertia  l39ibgAt2)'n  1.14(AgVt2)174 


Gravity-viscous 
Surface  tension-viscous 


1.39  (AgAV72  V172)'74 
1. 43(al  3p^v_1)1/4 


A  =  1  -  (po/  pw) 
v  =  v  of  water 

A  =  0.5  volume  of  oil  /  unit  length  of  oil  slick 
¥  =  volume  of  oil  slick. 


0.98(AgVV72  v-,/2),/6 

1  2,3  -2  -1)1/4 

1.60(af  pwv  ) 


a.  One-dimensional  spreading.  b.  Axisymmetric  spreading. 

Figure  8.  Definition  sketches  for  the  spreading  laws  for  oil  slicks. 


Fay's  spreading  theory  is  derived 
for  single-component,  constant-vol¬ 
ume  slicks  with  idealized  configu¬ 
rations  in  quiescent  water.  This  the¬ 
ory  considers  the  spreading  of  oil  to 
be  a  result  of  two  driving  forces — 
gravity  and  surface  tension — inter- 
balanced  by  the  retarding  forces  of 
inertia  and  viscosity.  The  spreading 
of  an  oil  slick  is  considered  to  pass 
through  three  phases.  In  the  begin¬ 
ning  phase,  only  gravity  and  inertia 
forces  are  important.  In  the  interme¬ 
diate  phase  the  gravity  and  viscous 
forces  dominate.  In  the  final  phase 
the  surface  tension  is  balanced  by 
viscous  forces.  Formulas  for  one-di¬ 
mensional  spreading  and  radial 
spreading  at  different  stages  are 
summarized  in  Table  2  using  terms 
defined  in  Figure  8. 

The  spreading  rate  during  each 
phase  can  be  obtained  by  taking  time 
derivatives  of  the  formulas  given  in 
Table  2.  The  equations  for  spread¬ 
ing  rates  are  summarized  in  Table  3. 


Table  3.  Spreading  rates  of  oil  slicks  and  phase  transition 
times. 

A.  Spreading  rates  for  one-dimensional  slicks,  dLJdt. 

Spreading  phase  Spreading  rate 

Gravity-inertia  (A g)'73  lL',/3(dL/df)  f273  +  L2/3  H/3] 

Gravity-viscous  139(AgV,/2),/4  l3/'  L  I-578  +  (dL/df)  f378! 

Surface  tension-viscous  1.43(op^  )172  (3/4  f  v)"174 

B.  Spreading  rates  for  circular  oil  slicks,  ARJAt. 

Spreading  phase  Spreading  rate 

Gravity-inertia  0.285fAg)' 74  (/' 72  V374)  [(dV/d/)  +  (2V/f )] 

Gravity-viscous  0.98fAgv-'  72),7‘  ('/,  V-273  f'74  (d¥/dt)  +  \¥'n  H74) 

Surface  tension-viscous  l^o(o2p^v"1C1)1/4 

C  Times  of  phase  transition 

Transition  One-dimensional  spreading  Radial  spreading 

Gravity-inertia  to  ( L 8  V3  A"2  g"2)’77  O-SSfVv’  A''  g-’)173 

gravity-viscous 

Gravity-viscous  to  (0.89 AgL4  o'2  pw2  v172)273  0.38(pw/a)(  AgvV2)1 73 

surface  tension-viscous _ 

V  «  volume  of  oil  per  unit  length  along  the  major  axis  of  the  slick  (L  is  a 
characteristic  length) 

¥  <*  total  volume  of  the  slick. 
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•  10,000  Barrels 
O  1,000 
■  100 


Time  (min) 

Figure  9.  Slick  radius  vs  time  (pw  =1.09  g/cm3,  po  =  0.9  g/cm3,  vw  =  2.543x 
10T2  cm1  Is,  a  =  11.02  dynes/cm).  The  points  represent  phase  transitions. 
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•  10,000  Barrels  - 
0  1,000 

■  100  - 

D  10 


luni 


innm 


Time  (min) 

Figure  10.  Slick  thickness  vs  time  (pa  =  1.09  g/cm3,  p0  =  0.9  g/cm3,  vw  = 
2.543xl(T2  cm2/$,  a  =11.02  dynes/cm).  The  points  represent  phase  transi- 


The  rate  of  change  of  the  oil  volume  in  these  equations  represents  the  change  due  to  weathering  and 
the  changes  in  oil  volume  distribution  in  various  parts  of  the  slick.  In  addition  to  the  rates  of  spreading 
at  different  phases,  the  times  at  which  each  phase  transition  occurs  also  need  to  be  determined.  These 
transition  times  can  be  obtained  by  letting  equations  in  the  appropriate  phases  equal  each  other  and 
solving  for  time.  Equations  for  the  transition  times  are  summarized  in  Table  3.  Fay  (1969, 1971) 


.6 


observed  that  the  changes  in  slick  properties  caused  by  weathering  may  result  in  the  eventual 
cessation  of  mechanical  spreading.  Based  on  a  number  of  field  observations,  Fay  suggested  that 

Af  =  105  V3/4  (32) 

where  A(  is  the  final  slick  area  in  m2  and  Y  is  the  total  volume  of  the  slick  in  m3.  In  this  study  the 
mechanical  spreading  is  considered  to  cease  when  the  slick  thickness  reduces  to  10-5  Y  1/4  m.  To  illus¬ 
trate  the  respective  regimes  of  the  different  phases  of  spreading,  variations  of  the  radius  and  thick¬ 
ness  of  circular  slicks  of  different  volumes  resulting  from  an  instantaneous  spill  are  presented  in  Fig¬ 
ures  9  and  10. 

The  formulas  in  Tables  2  and  3  were  derived  for  simple  slick  geometries  that  exist  under  idealized 
conditions.  In  the  simulation  model  the  radial  spreading  formulas  are  used  when  the  slick  is  nearly 
circular,  and  the  one-dimensional  formulas  are  used  when  the  slick  area  is  elongated.  A  slick  is  con¬ 
sidered  to  be  elongated  when  the  aspect  ratio  of  the  slick  area  is  greater  than  3.  The  aspect  ratio  refers 
to  the  length-to-width  proportion  of  the  slick,  and  the  orientation  refers  to  the  angle  0  between  the 
major  axis  of  the  slick  and  the  x  axis,  as  shown  in  Figure  1 1 . 


Figure  11.  Definition  sketch  for  variables  used  to  compute  slick  aspect  ratio  and  orientation. 

The  orientation  of  the  slick  is  computed  using  the  moments  and  product  of  inertia  of  the  slick.  The 
angle  0  can  be  expressed  using 

tan(20)  =  ^  (33) 

Ix-ly 

where 

/x=Zy2,  Jy=Ix2,  Pxy  =  Ixy  (34) 

where  I%  and  /y  are  the  moments  of  inertia  of  the  oil  slick  with  respect  to  the  x  and  y  axis,  respectively, 
and  is  the  product  of  inertia.  Once  the  orientation  is  known,  the  x  and  y  coordinates  of  all  the 
particles  in  the  slick  are  transformed  into  coordinates  in  an  x',  y7  coordinate  system  in  which  the  x' 
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axis  is  the  major  axis  of  the  slick.  The  aspect  ratio  is  computed  by  eq  35  based  on  the  x',  y7  coordinates 
of  all  the  oil  particles: 

5V 

Aspect  ratio  = - .  (35) 


This  aspect  ratio  is  used  to  determine  / 

whether  the  spreading  of  the  oil  slick  is  'Vx 

radial  or  one-dimensional.  The  use  of  .  ''/(f ' 

an  aspect  ratio  of  3.0  for  the  transition  /  /  1  “A 

from  radial  spreading  to  one-dimen-  ^  _  ,  I  Q 

sional  spreading  is  subjective,  but  it  ^  \  ^  H 

gives  reasonable  results.  /\  Centroid  of  _  — 

For  a  nearly  circular  slick  the  slick  h  \.  /  A  „ _ Slick 

area  is  divided  into  eight  pie-shaped  ""t  /  \ 

sectors  as  shown  in  Figure  12.  For  an  V,  /  O' 

elongated  slick  the  entire  slick  is  broken  1,11  \ 

up  into  a  series  of  strips  as  shown  in  Fig¬ 
ure  13.  For  both  cases  thespreadingrate 

ofeach  segment  is  calculated  using  Fay's  Figure  12.  Division  of  a  nearly  circular  slick  into  pie  segments. 
formula  independent  of  other  segments 

in  the  slick.  This  simplification  assumes  that  concentration  gradients  at  the  boundaries  between 
neighboring  segments  are  negligible  for  mechanical  spreading.  The  segmentation  also  allows  for 
different  spreading  rates  in  different  regions  of  the  slick,  providing  a  more  realistic  description  of  the 
field  situation.  During  each  time  step  the  increase  in  mean  radius  of  each  pie-shaped  sector  in  a  nearly 
circular  slick  or  the  mean  width  of  each  segment  in  an  elongated  slick  can  be  calculated  from  the 
spreading  formulas  in  Table  2. 


Figure  12.  Division  of  a  nearly  circular  slick  into  pie  segments. 


Strip 
Centroids 


Direction  of  Spread 


Figure  13.  Division  of  an  elongated  slick  into  strips. 


For  a  nearly  circular  slick  the  rate  of  outward  movement  of  an  oil  parcel  along  the  radial  direction 
in  a  particular  pie-shaped  sector  located  at  a  distance  r  from  the  centroid  of  the  slick  can  be  calculated 
from  the  spreading  rate  of  the  mean  radius  of  the  sector  as  (r/r)(dr  /dt).  Parcels  far  away  from  the 
main  slick  will  be  excluded  from  this  process  since  small  isolated  patches  of  oil  will  not  be  subjected 
to  mechanical  spreading.  In  this  study,  parcels  that  account  for  the  outer  5%  of  the  total  slick  volume 
are  excluded  from  the  mechanical  spreading  analysis.  This  is  equivalent  to  excluding  the  parcels 
located  at  a  radial  distance  greater  than  2 2r  from  the  centroid  of  the  slick.  This  was  determined  from 
numerical  experiments  for  the  spreading  of  circular  slicks  in  a  wide  rectangular  channel. 

For  an  elongated  slick  the  rate  of  outward  movement  in  the  width  direction  of  an  oil  parcel  in  a 
segment  with  mean  length  f,  located  a  distance  l  from  the  centroid  of  the  segment,  can  be  calculated 
from  the  rate  of  spreading  of  the  mean  length  of  the  segment  as  (1/  f  )(df  /dt). 
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Spreading  under  an  ice  cover 

The  spreading  of  oil  under  ice  covers  has  received  little  attention.  However,  the  study  of  Hoult  et 
al.  (1975)  has  provided  some  information.  Their  study  suggests  that  appreciable  mechanical  spread¬ 
ing  will  only  occur  during  continuous  spills.  The  oil  thickness  stabilizes  when  the  equilibrium 
thickness  for  the  flow  condition  is  reached.  There  are  no  pressure  gradients  or  surface  tension  forces 
to  cause  the  oil  to  spread  further.  The  oil  reaches  an  equilibrium  state  in  which  cavities  formed  by  the 
ice  roughness  contain  a  volume  of  oil  that  can  decrease  only  with  a  significant  increase  in  the  curren  t 
speed.  Since  a  continuous  spill  will  repeatedly  add  oil  to  fill  the  cavities,  the  excess  oil  will  effectively 
spread  to  the  empty  neighboring  cavities  and  establish  an  equilibrium  state  there.  This  is  a  crude  but 
reasonable  assessment,  since  only  the  excess  cil,  from  overfilling  the  cavities  under  the  ice  cover, 
would  be  expected  to  spread. 

Tire  formula  used  to  model  mechanical  spreading  for  a  continuous  spill  under  ice  is 

r  =  0.25  (^l|1/6f 2/3  (36) 

1  It'  I 


where  r 
A 

S 

Q 

h’ 


slick  radius 

(piv-po)/pw,  the  relative  density  ratio 
acceleration  due  to  gravity 

average  volume  flow  rate  since  the  beginning  of  the  spill 
half  of  the  root-mean-square  roughness  height  of  the  ice  cover. 


This  equation  is  a  result  of  balancing  the  frictional  drag  from  the  ice  cover  with  the  pressure  drop  that 
occurs  as  the  oil  flows  into  open  cavities.  In  the  simulation  models  no  mechanical  spreading  will 
occur  for  an  instantaneous  spill  or  once  the  continuous  oil  discharge  stops.  If  the  oil  discharge  is  in 
progress  and  the  slick  is  nearly  circular,  the  mechanical  spreading  will  be  calculated  by  eq  36. 

Shoreline  boundary  conditions 

An  oil  slick  will  reach  a  shoreline  sometime  after  a  spill  occurs.  Gundlach  and  Hayes  (1978)  pro¬ 
posed  a  method  for  classifying  shorelines  according  to  their  "vulnerability,"*  which  is  an  index  re¬ 
flecting  the  environmental  sensitivity  of  the  shoreline  to 
oil  pollution.  For  beaches  of  different  vulnerabilities  Tor- 
grimson  (1974)  suggested  the  use  of  "half-life"  values  to 
describe  the  ability  of  the  shore  to  retain  the  oil.  Half-life 
describes  the  "absorbancy"  of  the  shoreline  by  describing 
the  rate  of  re-entrainment  of  oil  after  it  has  landed.  Table 
4  presents  the  half-lives  for  different  types  of  shorelines 
along  with  their  vulnerability  indices. 

The  present  model  uses  the  half-life  concept.  The  vol¬ 
ume  of  oil  remaining  on  the  beach  is  related  to  its  original 
volume  by 

V2  =  ¥,  exp  [-k(t2-tj]  (37; 

where  ¥,  and  Y2  are  the  volumes  of  oil  on  the  beach  at 
times  fj  and  f2,  respectively,  and  k  is  a  decay  constant. 

Since  during  one  half-life  the  volume  of  the  oil  on  the 

*For  some  shorelines  in  the  United  States,  Environmental  Sensitivity  Index  maps  indicating  typesof  shoreline  characteristics 
are  available.  Maps  for  the  Detroit-St.  Gair  river  system  are  expected  in  the  near  future  (T.  Kaiser,  NOAA,  Ann  Arbor, 
Michigan).  This  information  can  be  used  with  the  current  model. 


Table  4.  Shoreline  descriptor  and  de¬ 
fault  parameters.  (After  Torgrimson 
1974.) 


Shoreline  Half  Vulnerability 

descriptor  life  index 


Exposed  headland 

1  hour 

1 

Wave-cut  platform 

1  hour 

2 

Pocket  beach 

1  day 

3 

Sand  beach 

1  day 

4 

Sand  and  gravel  beach 

1  day 

5 

Sand  and  cobble  beach 

1  year 

6 

Exposed  tide  flats 

1  hour 

7 

Sheltered  rock  shore 

1  year 

8 

Sheltered  tide  flat 

1  year 

9 

Sheltered  marsh 

1  year 

10 

Land 

1  year 

0 
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beach  will  be  reduced  by  half,  the  decay  constant  k  can  be  expressed  in  terms  of  the  half-life  X  as 


*  =  l-  ln(V2)]A.  (38) 

The  fraction  of  the  oil  re-entrained  into  the  water  body  during  each  time  step  is 

(V,  -  V2)/Vl  =  l-e~kAt=l-  0.5a'a.  (39) 

Evaporation 

Evaporation  can  account  for  the  largest  loss  in  oil  volume  during  the  early  stage  of  the  slick 
transformation.  In  this  study  the  formulation  developed  by  Mackay  et  al.  (1980)  is  used  to  calculate 
the  evaporation  rate  of  the  oil.  The  volume  fraction  of  oil  evaporated  is  determined  as 


F  =  (l/C)[lnPo  +  ln(CKEf  +  l/Po)]  (40) 

where  £  =  K£t  =  "evaporative  exposure"  term,  which  varies  with  time  and  environmental  condi¬ 
tions 

ke  =  KmAv/(RTVq) 

Km  =  mass  transfer  coefficient  (m/s),  0.0025^8d 
A  =  spill  area  (m2) 
v  =  molar  volume  (m3/mole) 

R  =  gas  constant  (82X10-6  atm  m3/(mol-K) 

T  =  surface  temperature  of  the  oil  (K),  which  is  generally 
close  to  the  ambient  air  temperature  T£ 

Vq  =  initial  spill  volume  (m3) 

Lfwlnd=  wind  speed  at  10  m  above  the  water  surface  (m/s). 

The  initial  vapor  pressure  Pq  in  atmospheres  at  the  temperature  TE  is 

lnPo  =  10.6[1  -(Tq/Te)]  (41) 

where  Tq  is  the  initial  boiling  point  (K)  and  T£  is  the  ambient  air  temperature  (K).  The  constant  C  can 
be  determined  by  the  relationship  T£C  =  constant.  C  values  for  T£  =  283  K  and  the  initial  boiling  point 
Tq  are  given  in  Table  5.  For  crude  oils  of  different  API  index  values,  C  values  at  T£  =  283  K  are  given 
in  Table  6,  along  with  T  .  This  table  can  be  replaced  by  the  following  functional  relationships  obtained 
through  curve  fitting  as  shown  in  Figures  14  and  15: 


Table  5.  Suggested  evaporation  para¬ 
meters  for  various  petroleum  fractions 
(Te  =  283  K).  (After  Mackay  et  al.  1980.) 


T  (K) 

O 

C 

P,  (atm) 

Motor  gasoline 
Summer 

314 

5.99 

0513 

Winter 

308 

6.23 

059 

Aviation  gasoline 

341 

2.81 

0.12 

Diesel  fuel 

496 

557 

3.4x10"* 

Jet  fuel 

418 

5.06 

6.0xl0-3 

No.  2  furnace  oil 

465 

7.88 

l.lxlO-3 

Lube  (heavy  and 

583 

8.61 

1.32X10"3 

light) 

Heavy  gas  oil 

633 

8.99 

2.0x10* 

Residuals 

783 

357 

755x10* 

Light  gas  oil 

473 

657 

8.1x10* 

Table  6.  Suggested  evaporation 
parameters  for  various  crude  oils 
(T.  =  283  K).  (After  Mackay  et  al. 
1980.) 


Specific 

gravity 


API 

(g/cm') 

C 

P  (atm) 

O 

10 

1.0 

89.2 

366 

0.044 

12 

0.986 

69.4 

348 

0.088 

15 

0.966 

52.1 

339 

0.13 

20 

0.934 

34.7 

329 

0.18 

25 

0.904 

27.2 

330 

0.17 

30 

0.876 

22.33 

325 

051 

35 

0.850 

195 

314 

051 

40 

0.825 

17.9 

304 

0.45 

45 

0802 

16.4 

283 

1.004 

20 


Coeff  icient 


from  Table  6. 


Figure  15.  Comparison  ofeq  43  with  T()  values 
from  Table  6. 


C  =  1158.9  API-1 1435 


(42) 


and 


T  =  542.6  -  30.275API  +  1 .565API2  -  34.39 API3  +  0.0002604API4 .  (43) 

O 

The  API  index  and  the  specific  gravity  of  the  oil  are  related  by 

Specific  gravity  =  141.5  /  (API  +  131.5) .  (44) 

The  molar  volume  of  oil,  which  is  required  in  eq  40,  can  vary  between  150X10"6  and  600x10"'’ 
m3/inole,  depending  on  the  composition  of  the  oil.  For  fuel  oils  this  value  is  approximately  200X10"6 
m3/inole.  The  molar  volume  can  be  computed  from  the  molecular  weight  of  the  oil.  Typical  molecular 
weights  for  various  oil  components  are  given  in  Table  7.  Tables  7  and  8  show  the  compositions  of 
several  different  oils,  along  with  the  molecular  weight  M(  of  each  component.  Based  on  these,  the 
molecular  volume  can  be  computed  from 

v  =  l /Kc(M()  (45) 

where  C.  is  the  weight  of  the  /th  component  per  unit  volume  of  oil. 
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Table  7.  Basic  data  for  oil  weathering.  (After  Moore  et  al.  1973.) 


Fraction 

Description  i 

%  by  wt 
in  crude  oil 

Density 

(g/mL) 

Boiling 

point 

CC) 

Molecular 

iveight 

Vapor 
pressure 
at  20 ‘C 
(mm) 

Solubility 

(g/rfg 

distilled  Hfi) 

1 

Paraffins 

«VcI2> 

0.1-20 

0.66-0.77 

69-230 

86-170 

110-0.1 

95-0.01 

2 

Paraffins 

<C.3-C25> 

0-10 

0.77-0.78 

230405 

184-352 

0.1 

0.01-0.004 

3 

Cycloparaffins 

«Vc12> 

5-30 

0.75-0.9 

70-230 

84-164 

100-1.0 

55-1.0 

4 

Cycloparaffins 

5-30 

0.9-1.0 

230405 

156-318 

1.0-0 

1.0-0 

5 

Mono-  and  di-cyclic 
aromatics 

(C6-C„) 

0-5 

0.88-1.1 

80-240 

78-143 

72-0.1 

1780-0 

6 

Polycyclic 

aromatics 

0-5 

1.1— 1.2 

240400 

128-234 

0.1-0 

12.5-0 

7 

Naphtheno- 

aromatics 

5-30 

0.97-1.2 

180400 

116-300 

1.0-0 

1.0-0 

8 

Residual 

(including 

heterocycles) 

10-70 

1. 0-1.1 

>400 

300-900 

0 

0 

Table  8.  Estimated  composition  and  comparison  of  solubilities  for 
various  petroleum  substances.  (After  Moore  et  al.  1973.) 


Composition  (%  by  xoeight) 


Fraction 

Description 

Crude 

A 

Crude 

B 

#2  fuel 
oil 

Kerosene 

Bunker 

C 

1 

Alkanes  (C6-C)2) 

1 

10 

15 

15 

0 

2 

Alkanes  (C^-Cj.) 

1 

7 

20 

20 

1 

3 

CycloparaffinsfCj-Cj  2) 

5 

15 

15 

20 

0 

4 

Cycloparaffins  (C13-Ca) 

5 

20 

15 

20 

1 

5 

Mono-  and  di-cyclic 
aromatics  (C6-Cn) 

2 

5 

15 

15 

0 

6 

Polycyclic  aromatics 

<c12-c18> 

6 

3 

5 

2 

1 

7 

Naphtheno- 
aromatics  (Cj-C^) 

15 

15 

15 

8 

1 

8 

Residual 

65 

25 

— 

— 

96 

Dissolution 

Dissolution  is  an  important  process  from  the  point  of  view  of  possible  biological  harm,  although 
it  accounts  for  a  negligible  fraction  of  the  mass  balance  of  the  oil.  Hydrocarbons  that  are  likely  to  dis¬ 
solve  in  the  water  are  also  likely  to  evaporate.  Dissolution,  which  tends  to  occur  in  the  first  few  hours 
of  a  spill,  thus  competes  with  evaporation.  Although  solubility  values  for  various  hydrocarbon  com¬ 
ponents  are  available,  these  values  are  difficult  to  use  correctly  for  modeling  purposes  since  they  are 
often  inconsistent  for  the  same  compounds. 


The  present  study  uses  the  method  of  Cohen  et  al.  (1980).  In  this  method  the  total  dissolution  rate 
N  is  calculated  by 

N  =  KAsS  (46) 

where  N  =  total  dissolution  rate  of  the  slick  (g/hr) 

K  =  dissolution  mass  transfer  coefficient,  assumed  to  be  0.01  m/hr 
=  slick  area  (m2) 

S  =  oil  solubility  in  water. 

Huang  and  Montastero  (1982)  suggested  that  the  solubility  for  a  typical  oil  can  be  calculated  by 
S  =  S  e~°  u  (47) 

O 

where  Sq  is  the  solubility  for  fresh  oil  and  t  is  the  time  in  hours.  Huang  and  Monastero  suggested  a 
typical  value  of  30  g/m3  for  Sq.  The  study  of  Lu  and  Polak  (1973)  provided  more  information  on  solu¬ 
bility;  they  formulated  the  rate  of  dissolution  as 

rd  =  cd  e-*  (48) 

where  rd  is  the  rate  of  dissolution  (mg  m2/day).  For  three  oil  samples  tested,  the  coefficients  c  and  d 
are  shown  in  Table  9. 


Table  9.  Dissolution  coefficients  at  25’C. 
(After  Lu  and  Polak  1973.) 


Oil 

type 

API 

c 

(mg  /nt2) 

d 

KS„=cd 
(gnr1  lir’) 

No.  2  fuel  oil 

35.5 

1043 

0.423 

0.0184 

Crude  oil 

38.6 

8915 

2380 

0384 

Bunker  Coil 

14.8 

459 

0303 

0.0104 

MODEL  APPLICATIONS 

Based  on  the  analytical  formulation,  two  computer  models  were  developed  and  applied  to  the 
connecting  channels  of  the  upper  Great  Lakes.  These  two  models  are  named  ROSS  (River  Oil  Spill 
Simulation)  and  LROSS  (Lake  and  River  Oil  Spill  Simulation).  ROSS  was  developed  for  simulating 
oil  spills  in  rivers  and  was  applied  to  the  Detroit  River,  St.  Clair  River,  lower  St.  Marys  River  and 
upper  St.  Marys  River.  LROSS  was  developed  for  simulating  oil  spills  in  lake-river  systems  and  was 
applied  to  the  Lake  St.  Clair-Detroit  River  system.  The  details  of  the  model  structure,  instructions  for 
input  data  preparation,  sample  output  and  program  listing  for  both  models  are  given  in  the  user's 
manuals  (Shen  et  al.,  in  prep.  a,b,  Yapa  et  al.,  in  prep.). 

To  illustrate  the  applicability  of  the  computer  models,  three  sample  simulations  are  presented 
here  (Fig.  16-18).  In  all  three  cases  no.  2  fuel  oil,  with  a  specific  gravity  of  0.84  and  a  surface  tension 
of  30  dynes/cm,  was  selected  as  the  spill  material.  The  other  input  parameters  for  each  case  are 
summarized  in  Table  10.  The  evaporation  characteristics  used  are  Tq  =  465  K,  C  =  7.88  and  v  =  2x1 0-4 
m3/mole.  The  solubility  of  the  oil  was  assumed  to  be  0.1873x1 02  Ib/ft3  A  time  step  of  15  minutes  was 
used  in  all  three  simulations. 
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Table  10.  Model  input  parameters  for  sample  simulations. 


Case  1 

Case  2 

Case  3 

Spill  location 

Lower  St.  Marys  River 

St.  Clair  River 
near  Frechette  Point 

Lake  St.  Clair 
near  St.  Clair 

Spill  type 

Continuous  leak, 
duration  =  30  min. 

Instantaneous 

spill 

Instantaneous 

spill 

Spill  volume 

35,000  L  (9,247  gal.) 

19,000  L  (5,020  gal.) 

37,850  L  (10,000  gal.) 

Oil  type 

No.  2  fuel  oil 
(sp.  gr.  =  0.84) 

Crude  oil 
(sp.  gr.  =  0.86) 

No.  2  fuel  oil 
(sp.  gr.  =  0.84) 

Row  conditions 

Discharge  =  111,017  ftVs; 
water  level  at  downstream 
side  of  Soo  Locks  =  581 .40ft 

Discharge  =  130,000  ftVs; 
water  level  at  upstream  side 
of  St.  Clair  River  =  575.00  ft 

Discharge  =  210,000  ftVs; 
water  level  at  lake 
outlet  =  575.04  ft 

Ice  conditions 

Open  water 

Partially  covered  downstream 
of  Detroit  Edison  with  15-cm 
ice  cover,  n,  =  0.035 

Open  water 

Wind  conditions 

0-3  hr  6.6  mph  from  W 

3-6  hr  6.6  mph  from  E 

6-14  hr  9.8  mph  from  N 

0-3.5  hr  1.3  mph  from  W 

0-24  hr  10  mph  from  NE 
24-40  hr  2  mph  from  E 

Diffusion  coefficient  (Du.) 

0.6 

0.6 

— 

Air  temperature  (°F) 

50 

20 

50 

Water  viscosity  (m2/s) 

1.308  x  KH 

1.792  x  10-6 

1.308  x  10-6 

Oil  viscosity  (poise) 

— 

0.158 

— 

Surface  tension  (dyne/cm) 

30 

30 

30 

Evaporation  parameters 

r„(K) 

465 

321 

465 

c 

7.88 

21.3 

7.88 

Solubility  parameters 

kSJ  g-2hr') 

0.0184 

0.884 

0.0184 

a  (day'1) 

0.423 

2380 

0.423 

SUMMARY  AND  CONCLUSIONS 

In  this  study,  two-dimensional  computer  models  for  simulating  oil  slick  movement  in  rivers  and 
lakes  were  developed.  The  models  were  then  applied  to  the  connecting  channels  of  the  upper  Great 
Lakes.  In  these  models  the  oil  slick  is  considered  to  be  a  collection  of  discrete  oil  patches.  The  transfor¬ 
mation  of  an  oil  slick  due  to  advection,  spreading,  evaporation  and  dissolution  are  considered.  In 
open-water  regions  the  advection  of  oil  patches  in  the  slick  are  determined  by  the  water  current  and 
wind  using  the  drifting  factor  formulation,  jn  ice-covered  regions  the  advection  of  the  oil  is  deter¬ 
mined  by  an  empirical  formula  developed  by  (T6xbh4  Schultz  (1981).  The  current  distribution  in  the 
lake  is  determined  by  a  rigid-lid  circulation  model  (Schwab  et  al.  1981).  The  current  distribution  in 
the  river  is  determined  by  a  stream-tube  method.  In  the  spreading  process,  mechanical  spreading  for¬ 
mulas  developed  by  Fay  (1969)  and  Hoult  (1972)  are  used.  TheseWmulas  consider  the  balance  of 
inertia,  gravity,  viscous  and  surface  tension  forces.  In  ice-covered  regions  the  formula  developed  by 
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Hoult  et  al.  (1975)  is  used.  In  addition  to  mechanical  spreading,  the  horizontal  turbulent  diffusion  of 
the  oil  patches  is  simulated  by  a  random-walk  formulation.  Formulations  developed  by  Mackay  et 
al.  (1980)  and  Cohen  et  al.  (1980)  are  used  to  determine  the  rate  of  evaporation  and  dissolution.  Boun¬ 
dary  conditions  along  the  shore  are  formulated  according  to  the  oil  retention  capability  of  the  shore¬ 
line. 

The  oil  slick  transformation  model  developed  in  this  study  contains  as  many  processes  <s  can  be 
effectively  and  analytically  modeled.  In  addition  to  computational  efficiency,  the  model  has  several 
special  features,  including  the  ability  to  model  instantaneous  and  continuous  spills,  the  ability  to  real¬ 
istically  describe  the  irregular  shapes  of  an  oil  slick  and  the  ability  to  account  for  the  time-dependent 
variation  of  the  flow  conditions.  Improvements  on  parts  of  the  model  will  be  possible  as  better  for¬ 
mulations  of  the  oil  slick  transformation  processes  become  available.  The  computer  programs  are  de¬ 
signed  so  that  it  will  be  easy  to  refine  the  model  elements  and  expand  the  model  to  include  additional 
slick  transformation  processes. 
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